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ABSTRACT 

The fundamental challenge in identification of nonlinear dynamic systems is determining the 
appropriate form of the model. A robust technique is presented in this paper which essentially 
eliminates this problem for many applications. 

The technique is based on the Minimum Model Error (MME) optimal estimation approach. 
A detailed literature review is included in which fundamental differences between the current 
approach and previous work is described. The most significant feature of the current work is the 
ability to identify nonlinear dynamic systems without prior assumptions regarding the form of the 
nonlinearities, in contrast to existing nonlinear identification approaches which usually require 
detailed assumptions of the nonlinearities. Model form is determined via statistical correlation 
of the MME optimal state estimates with the MME optimal model error estimates. The example 
illustrations indicate that the method is robust with respect to prior ignorance of the model, and 
with respect to measurement noise, measurement frequency, and measurement record length. 
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INTRODUCTION 


The widespread existence of nonlinear behavior in many dynamic systems is well- 
documented, e.g, Thompson and Stewart [1]; Nayfeh and Mook [2]. In particular, virtually 
every problem associated with orbit estimation, flight trajectory estimation, spacecraft dynamics, 
etc., is known to exhibit nonlinear behavior. Many excellent methods for analyzing nonlin- 
ear system models have been developed. However, a key practical link is often overlooked, 
namely: How does one obtain an accurate mathematical model for the dynamics of a particular 
complicated nonlinear system? General methods for actually obtaining accurate models for real 
physical systems are not nearly as widespread or well developed as are the techniques available 
for analyzing models. 

Accurate dynamic models are necessary for many tasks, including basic physical understand- 
ing, analysis, performance prediction, evaluation, life cycle estimation, control system design, 
etc. For example, most filter design assumes white process noise, yet many nonlinear effects are 
inherently non-zero mean; e.g., quadratic nonlinearities are always positive. In order to obtain 
a model with truly zero mean process noise for filter design purposes, all of the quadratic terms 
(and many other nonlinearities) must be well modeled. However, the complexity of many real 
systems greatly diminishes the possibility of accurately constructing a dynamic model purely 
from analysis using the laws of physics. 

Identification is the process of developing an accurate mathematical model for a system, 
given a set of output measurements and knowledge of the input. Many well developed and 
efficient identification algorithms already exist for linear systems (e.g., [3]-[7]). These often 
may be employed to model nonlinear systems when the system nonlinearities are small, and/or 
the system operates in a small linear regime. However, linearization does not work well (if 
at all) in every application, and even when it does provide a reasonable approximation, the 
approximation is normally limited to a small region about the operating point of linearization. 
Consequently, there is a real need for nonlinear identification algorithms. If nonlinearities are a 
predominant part of a system’s behavior, using a linear model to describe such a system leads 
to inconsistencies ranging from inaccurate numerical results to misrepresentation of the system’s 
qualitative behavior. Many important characteristics of nonlinear behavior, such as multiple 
steady-states, limit cycles, hysteresis, softening or hardening systems, chaos, etc., have no linear 
equivalent. Since nonlinearities are seldomly easily characterized, identification techniques may 
prove beneficial in developing accurate mathematical representations of nonlinear systems. 

Numerous methods for the identification of nonlinear systems have been developed in the 
past two decades. Many of these techniques are reviewed in Natke, Juang and Gawronski [8], 
Billings [9], and Bekey [10]. Most methods fall into one of the following categories: 

□ describing the nonlinear system using a linear model 

□ the direct equation approach 

□ representing the nonlinear system in a series expansion, and obtaining the respective coef- 
ficients either by using a regression estimation technique, by minimizing a cost functional, 
by using correlation techniques, or by some other approach 


444 


□ obtaining a graphical representation of the nonlinear term(s), then finding an analytical 

model for the nonlinearity 

With such diversity of nonlinear identification techniques, the choice of a particular algorithm 
may be based on criteria such as: the degree to which prior assumptions of the model form 
affect the user’s effort in applying the algorithms; the number of iterations required; the 
sensitivity to the presence of measurement noise in the data; the number of state measurements 
needed; whether or not knowledge of the initial conditions is required; the kind of forcing 
input(s) required or permitted (step, white gaussian noise, sinusoidal, etc.); the ability to handle 
hysteritic or discontinuous nonlinearities; the degree of a priori knowledge of system properties 
required; and the computational requirements. Most algorithms differ widely in at least some of 
these comparisons; the choice of a particular technique depends on the needs of the particular 
application. 

Among the methods which linearize the nonlinear system are those presented by Jedner and 
Unbehauen [11] and Ibanez [12]. Jedner and Unbehauen represent a nonlinear system, which 
may often operate in small regions around a number of operating points, by an equivalent number 
of linear submodels. It is assumed that the system operates at only a few points. Although the 
model may work well for controller design, the points at which the system is operating must be 
known and the linear models apply only within the operating regions. Ibanez takes a slightly 
different approach by assuming the system response to be periodic at the forcing frequency. An 
approximate transfer function is constructed. The transfer function is dependent on the amplitude 
as well as on the exciting frequency and is valid only within the region of exciting frequencies. 

The direct equation approach is used by Yasuda, Kawamura and Watanabe [13], [14]. The 
input and output measurements of a dynamic process are expressed in a Fourier Series using, for 
example, an FFT algorithm. The system nonlinearity is represented as a sum of polynomials with 
unknown coefficients. Applying the principle of harmonic balance, the polynomial coefficients 
as well as the other system parameters are obtained accurately. Knowledge of the nonlinearity 
is needed to construct the polynomial. Truncation in the Fourier Series expansion of the input 
or output may lead to error. 

The regression estimation approach is used by Billings and Voon [15] and Greblick and 
Pawlak [16]. Billings and Voon use the NARMAX model (Nonlinear Auto Regressive Moving 
Average model with exogenous inputs) to represent the nonlinear system. A stepwise regression 
method determines the significant terms in the NARMAX model. Then a prediction-error 
algorithm provides optimal estimates of the final model parameters. Greblick and Pawlak 
represent the linear dynamic submodel by an ARMA model and the nonlinearities by a Borel 
function. A non-parametric kernel regression estimation is employed to obtain the final analytical 
model. 

Kortman and Unbehauen [17] and Distefano and Rath [18] use the minimization of an error 
cost function as a means of obtaining the coefficients of the functions used to represent the 
nonlinearities. The method presented by Kortman and Unbehauen uses only system input and 
output information to estimate the polynomial representing the nonlinearities and the parameters 
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of the linear components. It is robust in the presence of noise, although iteration is necessary. 
Distefano and Rath present two techniques, a non-iterative direct identification and an iterative 
direct identification. In the first technique, measurement of all variables is required and the model 
parameters are obtained through the minimization of an error function. In the second technique, 
iteration is used to minimize a cost function yielding the system parameters in addition to the 
state trajectories. In Distefano and Rath, the nonlinear model form is also taken to be known. 

In other techniques, as in statistical linearization, a nonlinear relation is replaced by a linear 
equivalent gain. Broersen [19] extends the technique of statistical linearization by representing 
the nonlinearity as a linear combination of a number of arbitrary functions. Correlation techniques 
are then used to determine the coefficients of these functions. The number and type of functions 
selected depends on the desired accuracy as well as some knowledge of the system nonlinearity. 
Reasonable accuracy is obtained in the presence of noise and no iterations are necessary. 
Although some of the basic properties of the true nonlinear output are preserved, it is limited to 
only random excitation, and knowledge of all states and forcing terms is required. 

In the method of multiple scales (Hanagud, Mayyappa and Craig [20]), a perturbation solution 
to the nonlinear equation of motion is obtained. An objective function is built employing 
an integral least squares approach. The minimization of the functional yields the unknown 
parameters. Data on only one field variable is necessary, and the method is effective in the 
presence of high noise. The method of multiple scales, however, is restricted to systems 
with small damping and slight nonlinearities and, as in most other methods, the form of the 
nonlinearity is assumed a priori. The method typically requires some algebraic manipulations 
which may be quite involved, and these manipulations are only valid for a particular assumed 
nonlinear form. If the assumed nonlinear form is changed, the algebra must be repeated. 

Several techniques describe the nonlinear system using the Volterra or Wiener kernels. The 
Volterra series consists of the summation of impulse responses of increasing dimensionality. The 
Wiener series is also a set of orthogonal functions in which the input is white gaussian noise. 
Marmarelis and Udwadia [21], for example, estimate the first and higher order kernels appearing 
in the Volterra series using correlation techniques. Chen, Ishii and Suzumura [22] use cross- 
correlation functions in addition to the Volterra and Wiener series to describe nonlinear models 
and to show the relation between the system inner structure and the series. Although weakly 
nonlinear systems can be described by the first few kernels, for strongly nonlinear systems these 
series give accurate numerical results only at the expense of an excessive number of coefficients. 
This renders the analytical model impractical for control applications. 

Other popular series used in nonlinear identification are orthogonal polynomials such as 
Legendre (Wang and Chan [23]), Chebyshev, and Jacobi (Horn and Chou [24]). Horn and 
Chou expand the variables of the system into a shifted Jacobi series, reducing the nonlinear 
state equation into a linear algebraic matrix equation. The unknown parameters of the nonlinear 
system are then estimated using least squares. Even though the algorithm works well in the 
presence of noise, the nonlinear form must be known a priori. 
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Methods for the identification of nonlinear systems have also been developed based on 
the extended Kalman filter. The extended Kalman filter is the linear Kalman filter applied 
to nonlinear systems by linearizing the nonlinear model into a Taylor series expansion about 
the estimated state vector. Yun and Shinozuka [25] apply the extended Kalman filter for the 
parameter estimation of a quadratic term. The state vector is augmented by including the 
unknown parameters in addition to the state variables. Through a series of iterations, the response, 
as well as the unknown parameters, are estimated by the Kalman filter. Among its disavantages 
are high sensitivity to initial conditions, in particular if the initial conditions are barely known. 
The nonlinear form must be chosen a priori in order to estimate the corresponding parameter(s). 

Hammond, Lo and Seager-Smith [26] use an optimal control technique based on optimal 
control methods employed for linear system deconvolution. The form of the linear model is 
assumed to be known as well as the input and the output. A cost functional consisting of 
the weighted sum of the square of the error (between the actual and estimated output) yields 
an optimal estimated input. The estimated input and the actual input are used to obtain the 
nonlinearity as a function of the state variables. Although no previous assumption is made of 
the nonlinearities, there is no provision to deal with noise. 

All of the techniques outlined above have proven useful in certain applications. However, 
all of them are subject to one or more of the following shortcomings: 

1. The form of the nonlinearity (quadratic, cubic, exponential, etc.) must be assumed a priori. 
This is a very serious drawback, because the identification algorithm can only attempt to 
find the best model in the assumed form. If the form is assumed incorrectly, the resulting 
model may be so poor as to be useless, or it may appear to fit the data well enough that 
the user erroneously concludes that the correct model has been obtained. Also, for many 
techniques of this type, the effort required to test a given form is considerable, which greatly 
diminishes the effectiveness since multiple form tests are less likely to be conducted. 

2. Techniques which attempt to avoid the problem of a priori model form assumption through 
the use of series expansions generally eliminate any possibility of understanding the under- 
lying physics. Thus, although a good fit of the data might be achieved using a sufficient 
number of terms in the series, physical insight is lost. Moreover, large systems and/or par- 
ticularly complicated behavior may require that a very large number of terms be used to 
obtain a given level of accuracy. 

3. The presence of noise in the measurement data is not rigorously treated, yet noise is generally 
unavoidable. 

4. Initial conditions must be known in order to implement the algorithm. 

5. The algorithm can only be implemented if the data is obtained using very specific system 
excitations. 

The algorithm of the current paper compares favorably with existing algorithms in most of 
the categories listed above. It is robust with respect to measurement noise; does not require 
knowledge of the initial conditions; is independent of the forcing (but, like all methods, assumes 
that it is known); is not computationally prohibitive; and, most importantly, it requires minimal 
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a priori assumptions regarding the form of the model or the system properties. In fact, using the 
correlation technique outlined in the next section, the algorithm essentially eliminates the need 
to ever assume the nonlinear model form. 

The identification algorithm is based on a combination of Minimum Model Error (MME) 
state estimation, correlation techniques, and least squares. MME was first described by Mook 
and Junkins [27]. The MME combines the available measurements and an assumed model of 
the system to produce optimal estimates of the states and the model error. The assumed model 
represents an initial attempt to model the system using direct analysis, but may be extremely 
poor. Given the noisy output measurements of the system, MME estimates the state histories as 
well as the error in the assumed model. In previous work, the correct form and corresponding 
parameters of the nonlinear model were then estimated in a trial-and-error fashion, by assuming 
a nonlinear (in the states) form of the error terms, and then determining the best least-squares 
fit between the state estimates and the model error estimates. Thus, although the MME portion 
of the algorithm did not require the model form to be assumed, the subsequent least-squares fit 
between the state estimates and the model error estimates did. In Mook [28] it was shown that 
this approach could accurately identify terms in a Duffing oscillator, in the presence of noise and 
sparse measurements. The method worked well even when only a crude model of the dynamic 
system was assumed, and the error model used for the least-squares fit contained numerous terms 
in addition to the correct one(s). Later, in Mook and Stry [29], a simple harmonic oscillator 
with quadratic feedback was simulated on an analog computer. The algorithm was shown to 
accurately identify the nonlinear model from analog measurements. 

In this paper, the identification of the model from the MME-produced state and model error 
estimates is improved by using correlation techniques to select the form of the correction terms. 
The correction terms, when added to the initially assumed model, yield the true model of the 
system. The correction terms may consist of a combination of linear and nonlinear functions. 
An extensive library of linear and nonlinear functions has been assembled. The correlation 
technique is used to select the true forms from the library. Even when the true form of the 
nonlinearity was not present in the library, the correlation technique picks the closest form(s), 
typically, the first term(s) in the Taylor Series expansion. Once the forms have been selected by 
the correlation algorithm, least-squares is used to determine the model parameters. 

IDENTIFICATION ALGORITHM 

In this section, the identification algorithm is explained. First, the MME technique is briefly 
reviewed, and then the correlation technique used to automate the model form determination is 
explained in detail. 

The MME may be summarized as follows (a more detailed explanation may be found 
in Mook and Junkins [27]). Suppose there is a nonlinear system whose exact analytical 
representation is unknown, but for which output measurements are available. Using whatever 
means are available (analysis, finite elements, etc.), a system model is constructed. As shown 
in [27]-[29], the MME works well even if this system model is poor. The MME combines the 
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assumed model with the measurements to produce optimal estimates of (i) the state trajectories, 
and (ii) the error in the model. In the present work, these state and model error estimates are 
used for system identification. 

Consider a forced nonlinear dynamic system which may be modeled in state-space form 
by the equation 

x{t) = Ax(t ) + F(t) + /(x(f), i(<)) (!) 

where x(t) is the n x 1 state vector consisting of the system states, A is the n x n state matrix, 
F(t) is an n x 1 vector of known external excitation, and £(*(<), *(<)) is an n x 1 vector which 
includes all of the system nonlinearities. State-observable discrete time domain measurements 
are available for this system in the form 

]/(**) = £*(£(**)> 4) + to < tk < if ( 2 ) 

where y(tk) is an m x 1 measurement vector at time t^, is the accurate model of the 

measurement process, and y represents measurement noise, is assumed to be a zero-mean, 

gaussian distributed process of known covariance Rk- The measurement vector y(tk) ma y 
contain one or more of the system states. To implement MME, assume that a model, which 
is generally not the true system model because of the difficulties inherent in obtaining the true 
system model, is constructed in state-vector form as 

x(t) = Ax(t) + F(t) (3) 

Here, we show a linear model because in practice, linearization is the most common approach 
to modeling nonlinear systems. MME uses the assumed linear model in Eq. (3) and the noisy 
measurements in Eq. (2) to find optimal estimates of the states and of the model error. 

The model error, which includes the unknown nonlinear terms of the system, is represented 
by the addition of a term to the assumed linear model as 

x(t) = Ax(t) + F{t) + d{t) (4) 

where d(t) is the n x 1 model error to be estimated along with the states. 

A cost functional, J , that consists of the weighted integral square of the model error term plus 
the weighted sum square of the measurement-minus-estimated measurement residuals, is formed: 

pif 

+ / d(r) T W d{r)d.T (5) 

Jto 

where M is the number of measurement times, i{tk) is the estimated state vector and W is a 
weight matrix to be determined. 
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J is minimized with respect to the model error term, d{t). The necessary conditions for the 
minimization lead to the following two point boundary value problem (TPBVP), (see Geering 
[30]), 



(5u) 

A (() = -A t A(() 

m = -\wm 

(5b) 

(5c) 

H*t) = Hit) + mRt'Wk) - &(*&).<*)] 

( M) 

= !*(**).*» 


x(io) =X C or \(t 0 ) = o 

(Se) 

x(tf) = xj or A(f/) = 0 

( 5 /) 


where \(t ) is a vector of costates (Lagrange multipliers). Estimates of the states and of the model 
error are produced by the solution of this two-point boundary value problem. The estimates 
depend on the particular value of W. The solution is repeated until a value of W is obtained 
which produces state estimates which satisfy the “covariance constraint”, explained next. 

According to the covariance constraint, the measurement-minus-estimated measurement 
residual covariance matrix must match the measurement-minus-truth error covariance matrix. 
TTiis may be written as 


\^(h) - g k (Uh),tk)} T \i(t k ) - g k (Uh),tk)\ * R k ( 6 ) 

During the minimization, the weight W is varied until the state estimates satisfy the covariance 
constraint, i.e., the left hand side of Eq. (6) is approximately equal to the right hand side. The 
model error is, therefore, the minimum adjustment to the model required for the estimated states 
to predict the measurements with approximately the same covariance as the measurement error. 

The TPBVP represented by Eqs. (5a) to (5f) contains jumps in the costates and, consequently, 
in the model error. As evident from Eq. (5d), the size of the jump is directly proportional to 
the measurement residual at each measurement time. The noisier the measurements, the larger 
the jump size. A multiple shooting algorithm, developed by Mook and Lew [31], converts this 
jump-discontinuous TPBVP into a set of linear algebraic equations which may be solved using 
any linear equation solver. Multiple shooting also facilitates the analysis of a large number of 
measurements, by processing the solution at the end of every set of jumps. 

Correlation is a measure of the relationship that exists between two variables. The more 
highly correlated two variables are, the more closely will the change in one variable correspond to 
a change in the other variable. The cross-correlation coefficient between two discrete variables, 
say x and y, is defined as (see Newland [32] or Witte [33]) 


C(x,y) = Silfl — ~ y) 

CxO-yU 


( 7 ) 


450 


where n is the number of data points and the overbar denotes the mean of those n points. <r x 
is the standard deviation of the variable x and is defined as 



C( x ,y ) is a measure of the linear relationship between variables * and y. The value of C(x,t/) 
lies in the range -1 < C{x,y) < 1. If, for instance, changes in the value of * correspond to 
perfectly predictable (linearly) changes in the value of y, where the changes in both variables 
are of the same sign, then the value of C(x,y) is 1. If the changes are of opposite sign but 
still perfectly predictable, then the value of C(x,y) is -1. If changes in the values of x and y 
tend to correspond in sign but are not perfectly predictable, then 0 < C(x,y) < 1. If changes 
in the values of x and y tend to be of opposite sign but are not perfectly predictable, then 
-1 < C(x,y) < 0. If there is no linear relationship between the values of x and y, then 
C(x,y) = 0. For example, suppose x and y are multiples of each other, x = K *y, where K 
is an arbitrary constant of proportionality. Then 


C{x,y) 


S£i*K-f ll = 10 
EJ-i *0; - *) 2 ' 


(8) 


The true functional form of the model error can be found by calculating the correlation of the 
MME model error estimates with functions of the MME state estimates. If the functional form of 
the actual system is used, and if the estimates from MME are perfect, then C(x,y) = 1.0. Thus, 
an algorithm may be constructed which performs nonlinear system identification by (i) utilizing 
the MME to process the available measurements and the initial model in order to produce state 
estimates and model error estimates, and (ii) testing the correlation between the state estimates 
and the model error estimates usi'.g a “sufficient number” of functional forms so that the actual 
form is included among those tested. The MME does not require that the correct form of the 
model be known a priori. The correlation tests may be performed using an existing library of 
nonlinear functional forms, without input from the user. Thus, if the library is complete (in 
the sense that it contains the actual model form), the identification of the nonlinear model is 
accomplished, yet at no point in the algorithm is the user required to assume the correct model 

form. 

The success of the algorithm is determined by the ability of the MME to produce accurate 
state and model error estimates, and by the completeness of the library of nonlinear functions 
to be used in the correlation test. We now address these issues in order. 

The MME has been shown to consistently produce state and model error estimates of high 
accuracy in the presence of high measurement noise, low measurement frequency, and poor 
initial model [27-29]. Generally, however, some noise is still present in both the state estimate 
and the model error term, although these noise levels are considerably less than the noise in the 
original data. Let the model error term be given by x corTec tion = * + £ where f is the noise. 
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The cross-correlation between the error term and the test function y becomes 

EjLiQc* - g)(y» - y) + Ey=i £(yy - y) 


C (*,y) = 


1.0 


( 9 ) 


+ n £fc=l ( 2 £(* - *) + £ 2 ) 

As long as the noise is negligible all terms containing £ are small and affect the result only 
slightly. Thus, the correlation calculated for the actual function is close to, but not exactly equal 
to, 1, while the correlation calculated for incorrect terms remains close to 0. If the level of noise 
is excessive, say, of comparable magnitude to one or more of the actual nonlinear model terms, 
then the ability of the correlation test to distinguish this term from similar terms may be greatly 
reduced or eliminated. However, subsequent least-squares fit of the terms has, in every case 
tested, correctly selected the actual nonlinear function from among those which the correlation 
test could not distinguish. An example of this is shown in the next section. 

The issue of completeness of the library is now addressed. The error term may be composed 
of more than one function from the library, or the actual function may be missing from the 
library. Consider first the case where the actual error is a combination of library terms, say, two 
terms. The error term may be written Xco TTec ti on = xj -f-2 2 and the cross -correlation has the form 

E,-=i(sit ~ fiKw ~ y) + ]Cy=i( x 2 j — ^2 )(yy — y) 


C(x,y) = 


( 10 ) 


+ 0-z, 2 + I ELl 2 (*lJb + *l)(®2* + *2) 

The cross-correlation is highest for the term which constitutes the largest part of the error. Thus, 
it is desirable to execute the algorithm iteratively. The library term which constitutes the largest 
portion of the actual model error is identified first and then added to the MME model. The 
entire process (including MME) is then repeated, so that new state and model error estimates are 
obtained (note that the change in state estimates should be minimal, while the change in model 
error estimates should be a large reduction in magnitude). The largest term remaining in the 
model error is identified in each pass, then added to the initial MME model. 

An alternative to iterative application of the algorithm is to test the correlation of combi- 
nations of the library functions. An algorithm can be constructed which tests every possible 
combination of the functions explicitly contained in the library. This approach has not been 
attempted in the examples which follow. 

If the actual model error is not present in the library, then test cases show that the highest 
correlation values are calculated for the terms in the series expansion of the actual function. 
Thus, for example, if the actual model error was of the form jtn(x), but sin(x') was not present 
in the library, the correlation coefficients are highest for the terms x, x 3 , x 5 , etc. However, the 
test described by Eq. 7 is very fast, so the library may contain a very large number of terms. 

The final step in the identification procedure is to use a least-squares algorithm to fit the 
model error to the functional forms (i.e., perform parameter identification once the true nonlinear 
form has been determined). The error term is expanded into a combination of the functional 
forms such as 


d(t) = afx(x(t)) + 0f 2 (x(t)) + 7 /s(*( 0 ) + • • • 


( 11 ) 
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where a, (3, 7, ... are unknown coefficients to be determined by least squares, and /1, fi, 
f$, . . . are functions which are selected as a result of the correlation test (often, however, only 
one function is used at a time). Other parameters may be present inside the functions (such 
as, for example, coefficients of exponents). Eq. (11) may be sampled repeatedly (using the 
MME estimates) to obtain 

d(t 1 ) = a/i(x(fi)) + /?/ 2 (x(fi)) + 7/3(*(fi)) + • ■ • 

d(t 2 ) = a/i(x(t 2 )) + /?/2(®(*2)) + 7h{*{h)) + ■ • • 

d(ti) = afi(x{t[)) + Pf 2 (x{ti)) + + ... 

or, in matrix form, 

fix 1 = M lxr P pxX (12) 

where P = [a /3 7 . . ,] T is the vector of coefficients for the terms in d(t). Since estimates of 

d(t ) are available continuously throughout the time domain, the parameter l may be chosen quite 
large to improve the least squares fit. Generally, because of the potential jump discontinuities 
in the model error estimates at the measurement times, it is desirable to pick the least squares 
sampling times in Eq. (12) at points other than the measurement times. The least squares 
estimate is found by minimizing the following cost functional with respect to P. 

$ = [D - MP} t [D - MP] (13) 

The solution is given by 

p = (M t M)~ 1 M t R (14) 

If the functions include parameters to be estimated, the equivalent nonlinear least-squares problem 
is constructed. 

The multiple shooting algorithm presented by Mook and Lew [31] was used to obtain the 
MME solutions used in the tests presented in this paper. It was assumed in the examples that 
MME obtained the dynamic error term without knowledge of the boundary conditions on x, 
so some distortion of the correction term at the initial and final times was expected due to the 
constraints of Eqs. (5e-5f), i.e., by assuming no state knowledge is available at <o or if, we 
constrain A(< 0 ) = 0 and A (</) = 0. Therefore, in all test cases, the initial and final ten percent 
of the correction term data was ignored in the least squares fit. 

EXAMPLES 

For illustrative purposes, the true system was chosen as a simple harmonic oscillator with 
various forms of nonlinear feedback. The true system can be modeled as 

(i) = (-°i 0 ) (*) + (/(“«)) (15) 
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where x is position, v is velocity and the dot indicates differentiation with respect to time. For 
simplicity, the system was unforced. The term f(x,v) represents the nonlinear terms to be 
identified by the MME-based identification algorithm. Measurements were generated from the 
true system, Eq. (15), with different kinds of nonlinear functions f(x,v). The ability of the 
identification algorithm to identify the model with no prior knowledge of f(x,v'j is tested. Table 
1 shows the functions used in each simulation. Note that the unknown error term may be a 
combination of linear and nonlinear functions. Table 1 also shows the initial conditions and the 
amount of noise used to generate measurements for each test. The noise levels represent the 
percentage of the peak system response (actual percentages are higher for the majority of the 
measurements since the response is only at peak amplitude for brief periods). 

Table 1 

SUMMARY OF TEST CASES 


TEST# 

TRUE ERROR: f(x,v) | x(0) 

v(0) 

NOISE 

1 

3.0*x*x 

0.175 

0 

0 

2 

-0.1*x*x*v 

0.175 

0 

0 

3 

-0.5*cos(x)*cos(v) 

0.175 

0 

0 

4 

-1.0*v*sin(x) 

0.175 

0 

0 

5 

-1.0*x*x - 0.25*v 

0.350 

0 

0 

6 

-1.0*x*x*x - 0.1*tan(v) 

0.873 

0 

0 

7 

-1.0/cos(x) - 1.0*sin(v) 

1.750 

0 

0 

8 

3.0*x*x 

0.175 

0 

10% 

9 

-1.0*x*x - 0.25*v 

0.350 

0 

10% 

10 

-1.0*x*x*x - 0.1*tan(v) 

0.873 

0 

10% 


The assumed model used for the MME analysis consisted of the undamped linear oscillator 
part of the system, 

(:)-(-.;)(:) <■•> 
For each test, 200 measurements of position were obtained from the digital simulation of Eq. (15) 
at a sampling rate of 10 Hz. The functional form of the dynamic error, f(x,v), was determined 
solely from the least-squares fit of the functions identified during the correlation tests on the 
MME state and model error estimates obtained using only the model in Eq. (16). 

A library of functions was built consisting of approximately 300 of the most commonly 
found nonlinear and linear forms. For a particular test, after the model error term was found 
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from MME, it was correlated with each one of the functions in the library. The correlation test 
of the entire library of functions did not take more than a few seconds to execute, since the 
calculations are simple. The functional form of the unknown nonlinear term was chosen as the 
one for which the absolute value of the cross-correlation coefficient was closest to 1. Table 2 
shows the results for all 10 tests, including the true dynamic error, the highest cross-correlation 
coefficient obtained, the corresponding functional form, and the respective coefficient computed 
from the least squares fit. The star (*) indicates tests performed from noisy measurements. 

Table 2. 

IDENTIFICATION RESULTS FOR EACH TEST CASE 


TEST# 

TRUE ERROR(S) 1 C(d(t),fTI 

SELECTED 

L.S. 

1 

3.0*x*x 

0.999 

x*x 

2.99 

2 

-0.1*x*x*v 

0.999 

x*x*v 

-0.10 

3 

-0.5*cos(x)*cos(v) 

0.999 

cos(x)*cos(v) 

-0.49 

4 

-1.0*v*sin(x) 

0.999 

v*sin(x) 

-1.00 

5 

-1.0*x*x 

0.999 

x*x 

-0.99 

-0.25*v 

0.746 

V 

-0.24 

6 

-1.0*x*x*x 

0.936 1 

x*x*x 

-1.00 

-0.1*tan(v) 

0.999 

tan(v) 

-0.10 

7 

-1. 0/cos (x) 

0.927 

l/cos(x) 

-0.99 

-1.0*sin(v) 

0.999 

sin(v) 

-1.00 

8* 

3.0*x*x 

0.797 

x*x 

3.12 

9* 

-1.0*x*x 

0.937 

x*x 

-0.90 

-0.25*v 

0.772 

V 

-0.22 

10* 

-1.0*x*x*x 

0.838 

x*x*x 

-0.98 

-0.1*tan(v) 

0.583 

tan(v) 

-0.10 


For tests 1, 2, 3, and 4, the exact form of the nonlinearity was contained in the library and 
the measurements did not contain noise. The calculated value of C(d(t),f) was 1 for the true 
forms. In test 8, the library contained the exact form of the nonlinearity but the measurements 
contained significant noise. The correlation for the correct term was much higher than for any 
other term, but was approximately 0.8 instead of 1 due to the noise. In the cases where the 
error term consisted of two functions but the measurements were noise-free (tests 5, 6 and 
7), C(d(t),f ) was close to one for both functions after applying the algorithm iteratively as 
described in the previous section. 

When noise and more than one function was present in the dynamic error term (tests 9 and 
10), the maximum value of the cross-correlation coefficients dropped significantly and in some 
cases did not immediately identify the actual form over other similar forms. As an example, 
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Table 3 shows the top five cross-correlation values for the identification of the tan(v) term in 
test case 10. Note that the functions with the highest cross-correlation values are all similar in 
form to tan(v), and the corresponding correlation coefficients are of similar magnitude. Since 
^(^(0*/) not clearly identify tan(v ) as the missing term, the five functions yielding the 
highest C(d(t), f) values were individually least-squares fit to the model error term. In all cases 
(i.e., repeating this test for a number of different random noise samples), the function with the 
smallest least squares error cost was the correct function (fan(v)). Thus, the least-squares fit 
of the parameters to the functional forms also serves as a second test if the correlation test is 
inconclusive due to high noise levels. 


Table 3. 

HIGHEST CROSS-CORRELATION COEFFICIENTS 
OBTAINED FOR THE TAN(V) TERM OF TEST CASE 10 


FUNCTION 

C(d(t),f) 

L.S. 

L.S. cost 

tan(v) 

0.583 

-0.104 

0.588 

V 

0.584 

-0.119 

0.623 

v*cos(x)*cos(v) 

0.584 

-0.150 

0.659 

v*cos(x) 

0.586 

-0.126 

0.607 

sin(v)*cos(x) 

0.586 

-0.133 

0.621 


The number of data points used in the MME algorithm was irrelevant as long as there were 
enough points to reasonably span the qualitative aspects of the system (e.g., sinusoidal terms 
cannot be identified if the data only spans a small fraction of the period). 

If the exact functional form of the dynamic error term was not in the function library, the 
correlation procedure would pick the first term in the Taylor Series expansion of the exact form. 
For example in a test case where the dynamic error term corresponded to x*sin(v ) and x*sin(v) 
was deleted from the library, the function with the largest C(d(t),f) was x*v. Similarly, in 
several examples which are not shown the magnitude of the states, x and v, were small. Thus, 
the trigonometric functions of position and velocity were approximately equal to the first term in 
their Taylor Series expansions, i.e., coa(:c) « 1.0, st'n(a:) « x, cos(v) « 1.0 and sin(v ) « v. In 
these cases, assumptions of linearity are clearly valid, and are not of interest in the present work. 


SUMMARY AND CONCLUSIONS 

In this paper, an algorithm based on the MME estimation technique, coupled with correlation 
tests and least squares, has been developed for identification of nonlinear systems. The results 
of the examples indicate that the correlation technique applied to the MME-produced state and 
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model error estimates enables the form of the model to be accurately determined, thus eliminating 
the requirement that the form be assumed a priori. Once the form is determined, the least-squares 
fit provides excellent parameter identification. In cases of high noise, where the correlation test 
may not be able to distinguish the actual form from similar forms, the least-squares fit also 
proved to be a reliable second test for determining the actual form. 

At no point in the algorithm is the user required to assume the form of the model, representing 
a tremendous advantage over existing techniques, including the previous MME-based work. The 
MME does not require an accurate model in order to produce accurate state and model error 
estimates, and the correlation tests are automatically performed on a large existing library of 
functions. Additional functions and more sophisticated methods of combining existing functions 
can be added to the correlation testing portion of the algorithm (the authors are currently pursuing 
this), virtually eliminating the likelihood that the actual model error terms are not tested. 
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